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ABSTRACT 

The macroscopic behaviour of cosmic rays in turbulent magnetic fields is discussed. 
An implementation of anisotropic diffusion of cosmic rays with respect to the magnetic 
field in a non-conservative, high-order, finite-difference magnetohydrodynamic code is 
discussed. It is shown that the standard implementation fails near singular X-points 
of the magnetic field, which are common if the field is random. A modification to 
the diffusion model for cosmic rays is described and the resulting telegraph equation 
(implemented by solving a dynamic equation for the diffusive flux of cosmic rays) is 
used; it is argued that this modification may better describe the physics of cosmic 
ray diffusion. The present model reproduces several processes important for the prop- 
agation and local confinement of cosmic rays, including spreading perpendicular to 
the local large-scale magnetic field, controlled by the random-to-total magnetic field 
ratio, and the balance between cosmic ray pressure and magnetic tension. Cosmic ray 
diffusion is discussed in the context of a random magnetic field produced by turbulent 
dynamo action. It is argued that energy equipartition between cosmic rays and other 
constituents of the interstellar medium do not necessarily imply that cosmic rays play 
a significant role in the balance of forces. 



1 INTRODUCTION 

' The importance of cosmic rays for the dynamics of the in- 
terstellar medium (ISM) has long been recognized dParkerl 
! ll966l:lBerezinskii et al]ll99d) . Spatial gradients of the cos- 
' mic ray pressure contribute significantly to the force balance 
] in the ISM. If cosmic rays are confined within magnetic flux 
■ tubes, then the tendency toward pressure equilibrium re- 
' duces gas pressure within the tubes. Depending on the ef- 
, ficiency of cooling, either temperature or entropy will be 
approximately uniform across the tube, but in both cases 
density inside the tube will be decreased relative to the ex- 
terior, making the tube buoyant. This process is similar to 
magnetic buoyancy. Therefore, cosmic rays facilitate disc- 
halo connections in spiral galaxies by enhancing the buoy- 
ancy of magnetic structures in the interstellar gas. In the 
sun, magnetic buoyancy drives magnetic fiux tubes to the 
surface to form bipolar regions. In galaxies, magnetic buoy- 
ancy is believed to be strongly assisted by cosmic rays. 

The effects of cosmic-ray driven buoyancy are believed 
to be important for the opera tion of the galactic dynamo 
JParkeilll99l iMoss et al.lll999t) . This can help to speed up 
the growth of the magnetic field and maintain strong field 
amplification and regeneration, especially in the non-linear 
regime (Hanasz ct al. 2004). Many studies of the Parker in- 
stability as well as recent simulations of the galactic dy- 
namo rely on a hydrodyna mic description of cosmic rays 
JSchlickeiser fc LerchdlTgSSI) . which is especially convenient 
in models involving the large-scale dynamics of the interstel- 
lar medium. In this approach, the cosmic ray transport equa- 
tion for the phase space distribution function is integrated 



over particle momenta which results in a hydrodynamic-type 
equation for the cosmic ray energy density or pressure. Our 
aim here is to use this approach in order to clarify the re- 
lation between cosmic ray energy density and properties of 
the interstellar medium. 

Energy equipartition (or pressure balance) between cos- 
mic rays and magnetic fields is a common assumption used in 
radio astronomy, where it is used to estimate magnetic field 
strength from synchrotron intensity. A physical basis for this 
idea remains elusive and only qualitative arguments, related 
to cosmic ray confinement by magnetic fields, are used to jus- 
tify this concept. The assumption comes into question since 
the spatial distribution of cosmic rays may not precisely fol- 
low that of magnetic field strength. Furthermore, the idea 
of overall (statistical) pressure balance in the ISM would 
be more difficult to maintain if both magnetic and cosmic 
ray pressures are enhanced or reduced at the same posi- 
tions simultaneously. Recent arguments of iPadoan fc Scalol 
1:2005) suggest that, if the streamin g velocity of c osmic raw 
is proportional to the Alfv en speed jPelice fc Kulsrud,.200ll : 
iFarmer fc GoldreichI l20o3. and references therein), the lo- 
cal cosmic ray density is independent of the local magnetic 
field strength, and rather scales with the square root of 
the (ionized) gas density. Indeed, if both the magnetic fiux 
and the cosmic ray flux are conserved, BS = const and 
ricUS = const (where B is the magnetic field strength, S 
is the area within a fiuid contour, ric is the number density 
of cosmic rays and U is their streaming velocity), one ob- 
tains ricU/B = const, which yields Uc oc nV^, given that 
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[/ = Va OC BtI; , with rn the ion number density and Va 
the Alfven speed. 

We use a two-fluid model, where cosmic rays are de- 
scribed by an equation for their pressure (or energy density) 
and an equation of state. The cosmic rays are assumed to act 
directly on the background gas via their pressure gradient. 
We do not include any explicit means of exciting hydromag- 
netic waves by cosmic rays leading to t heir confinemen t (for 
a discussion of confinement issues, see ICesarskvllToSol) . but 
instead parameterize these processes by choosing an appro- 
priate advection velocity (as a superposition of the gas and 
Alfven velocities). There are several interesting questions re- 
garding high energy cosmic rays and their acceleration (e.g. 
Hillas 2005), which we are not attempting to address here. 
Instead, we want to know which process is mainly respon- 
sible for limiting the cosmic rays energy density and what 
is the relation of cosmic ray energy density with the mag- 
netic field. Is there local equipartition, or is there only global 
equipartition on the scale of the galaxy? Finally, we are in- 
terested in studying those effects in the ISM dynamics that 
only arise in the presence of cosmic rays. We begin with 
the governing equations and discuss issues that arise in con- 
nection with the numerical implementation of cosmic ray 
diffusion along magnetic field lines. 



where Kij is the diffusion tensor. The latter can be written 

as 

K,j = K^5,j + (if|| - K^)B,Bj, (8) 

where B = B/\B\ is the field-aligne d unit vector (e.g. 
iBerezinskii et alllQQCtlHanasz fc Lesch 20031. Here, K^^ and 
K± are the cosmic ray diffusion coefficients along and per- 
pendicular to the field, respectively. 

We assume ideal-gas equations of state for both the cos- 
mic rays and the gas, i.e. pc = (7c — l)ec andpg = (7g — l)eg, 
where 7c and 7g are the ratios of the total number of degrees 
of freedom to the number of translational degrees of freedom 
for the cosmic rays and the gas. Unless stated otherwise, we 
assume 7c = 4/3 a nd 7g = 5/3. Oth er choices for 7c include 
5/3 and 14/9 fe.g. IRvu et aljl2003l . and references therein). 

The system can be driven by an external force / in the 
momentum equation ||1J, and F in that equation includes 
additional forces such as the viscous force, V • {2i'pS), where 
u is the viscosity and Sij = ^(itij -I- — \5ijUk,k is the 
traceless rate of strain tensor, where commas denote partial 
differentiation. Furthermore, Qk ~ 2pvS'^ and Qm = rj/ioJ^ 
denote the viscous and Joule heating, and Qc is a cosmic 
ray energy source. 



2 METHOD 

2.1 Basic equations 

The hydromagnetic equations, supplemented by the 
advection- diffusion equation for the cosmic ray energy den- 
sity, and the cosmic ray pressure in the momentum equation. 



are 

1^ + V ■ ipu) = 0, (1) 

% + (ccw)+PcV- It = Dc + Qc, (2) 

^-HV- (egu)+j3gV-u = Dg + Qk + Qm, (3) 

^ + V-{puu) + V{p^+Pc)=JxB + f + F, (4) 

— = V X (u X B - ntJ-oJ), (5) 



where p, u and Pg are the gas density, velocity and pressure; 
Cc and Pc are the cosmic ray energy density and pressure, B 
is the magnetic field, J — 'V x B/pt} is the electric current 
density, rj is the magnetic diffusivity, Dg — V ■ (K'VT) is the 
thermal diffusion term (treated here isotropically; thermal 
diffusion is unimportant in the present context, but weak 
diffusion is necessary for numerical reasons). Further, T is 
the temperature related to the internal energy density (per 
unit volume), eg, via eg = pCvT, and Dc is the divergence of 
the diffusive cosmic ray energy flux taken with the opposite 
sign, i.e. 

Dc = -V- Tc. (6) 

The usual approach is to treat this term as Fickian diffusion, 
i.e., to assume that the flux is proportional to the instanta- 
neous gradient of the cosmic ray energy density, 

Tci ~ —KijVjEc (Fickian diffusion), (7) 



2.2 Non-Fickian difFusion 

Typical values of the diffusivit y along the magnetic fi eld are 
of the order 10^* cm^ s"^ Ce.g. lBerezinskn et alj|l990l) . Such 
large values would severely limit numerical modelling since a 
large diffusivity requires that the computational time step is 
small to ensure numerical stability; for example simulations 
with a resolution of 1 pc would require a time step of 10 years 
or less (e.g., Hanasz & Lesch 2003 reduce K^^ by a factor of 
10 to make the system tractable numerically) . This problem 
could be circumvented by employing an implicit numerical 
scheme. In the context of cosmic ray propagation, one would 
expect the advection speed to be not too much larger than 
the Alfven speed. Before discussing a possible remedy to this 
problem we note that, in the case of field-aligned diffusion, 
the problem can be even more severe. If we use the product 
rule and write Dc — Vi{Kij\7 jBc) in the form 

Dc = -Uc■'Vec + K^J^^^Jec, (9) 

we see that Ud ~ —dKij /dxj plays the role of a velocity 
transporting cosmic rays perpendicular to curved field lines. 
This term is proportional to the divergence of the dyadic 
product of unit vectors, V • (BB). At magnetic X-points, 
this term is singular, as explained below (we note that 0- 
type singular magnetic points do not cause difficulties). 

We illustrate this complication using a simple magnetic 
field configuration B = (a;, — t/, 0)"^ with a null point at the 
origin, which leads to the singular behaviour of V • {BB), 
and hence to a singularity of \Uc\- 

V ■ (BB) = -J ( i3x' - y')y j , (10) 

where = x'^ + y"^ . This expression diverges at the origin 
and leads to infinite propagation speed which would, tech- 
nically speaking, limit to zero the length of the time step of 
an explicit timestepping scheme. In spite of this singularity. 
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the cosmic ray energy density must stay finite. In fact, one 
can sliow that, in a closed or periodic domain, the maximum 
cosmic ray energy density, max(ec), can only decrease with 
time. This is a well-known general property of the diffusion 
operator; in Appendix^ we derive this result for the form 
of the diffusion tensor appropriate for cosmic rays. The rea- 
son that majc(ec) can remain finite, despite V • (BB), and 
hence Uc, becoming infinite, is that the parabolic system 
of equations can adjust itself instantaneously so that Vcc 
tends to zero where Uc diverges. A numerically convenient 
remedy to this problem will be discussed in Sect. 12.31 where 
a non-Fickian diffusion model is used. In the following we 
describe this approach in more detail. 

A physically appealing, and widely adopted way to im- 
prove the diffusion equation so as to limit the propagation 
speed to a finite value involves a more accurate descrip- 
tion of the diffusive flux. This generalization has been ap- 
plied, e.g., to turbulent diffusion. In turbulence, the clas- 
sical turbulent diffusion equation, dn/dt = Dd'^n/dx^, 
arises if the turbulent velocity field is assumed to be 5- 
correlated in time; this approximation is consistent with 
Eq. (|7J or its simplifications. In order to ensure finite prop- 
agation speed of the diffusing substance, it is sufficient to 
allow for a finite correlation time r of the velocity field. 
This leads to equation (IIH for the diffusive fiux. The cor- 
responding equation for the diffusing quantity reduces to 
the telegraph equation dn/dt + rd^n/dt^ — D d^n/ dx^ , or 
its generalizations. These arguments have been recently dis- 
cussed by iBakuninI (|2p03a b) . The telegraph equation has 
been used to correct acau sal cosmic-ray diffusion models 
(e.g., iGombosi et"ai]|l993h . This type of non-Fickian dif- 
fusion also eme rges quite na turally in turbulent diffusion of 
passive scalars jBlackman fc Field 2003) and has been con- 
firmed in direct simulations (Brandenbur g et al.ll20oJ^ . On 
long enough time scales, or for sufficiently small values of 
r, the non-Fickian description of diffusion reduces to the 
Fickian limit. 

Thus, we replace Eq. O by 



T 



(non-Fickian diffusion), (11) 



dt 

where Kij — rKij corresponds to the original diffusion ten- 
sor. Similarly to Eq. (|5J, we write 



(12) 



Quantitatively, the deviation from Fick's law is controlled 
by the dimensionless parameter 



St = 



K 



1/2. 



(■^11 



Nl/2 



(13) 



where £ is the typical length scale of the initial struc- 
ture. In the context of turbulent diffusion, this dimension- 
less parameter is often referred to as the Strouhal num- 
ber jLandau fc LifshitdlToSTl: iKrause fc Radlerlll980t) . The 
larger the Strouhal number, the more important are non- 
Fickian effects resulting in a wave-like behaviour of the so- 
lution. Unlike the solution of the classical diffusion equation, 
where an initial perturbation to the trivial solution has an 
effect at every position for any t > 0, solutions with non- 
Fickian diffusion remain unperturbed ahead of a propagat- 
ing front. 
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Figure 1. The spread of an initial Gaussian distribution of cos- 
mic ray energy density (of a half- width £): the distribution at a 
time t = 1 is shown, as a function of x/£, for three values of 
the Strouhal number St. Note that the behaviour of the solution 
becomes more wave-like as St increases. 



A suitable estimate of the Strouhal number can be ob- 
tained assuming that the relevant correlation time is of the 
order of the Alfven crossing time for magnetic structures of 
scale £, i.e. St ~ {K^^/Va£)^^'^- This yields (for gas number 



density 0.1 cm ) 



St ~ 20 



B 
5nG 



10 pc 



(14) 



^4 X 1028 cmVs 

In Fig. we illustrate the one-dimensional spread of 
an initial Gaussian distribution of cosmic rays, ec = 
exp{ — ^x^/£^) after t = t for three values of St. For small 
values of St, the solution evolves similarly to that of the dif- 
fusion equation (solid and dotted lines in Fig. . For large 
values of St, the distribution of cosmic rays develops two 
local maxima of ec that propagate outwards as shown with 
dashed line, a typical wave-like behaviour. In the limiting 
case of very large values of St the governing equation re- 
duces to the wave equation, and the classical diffusion is 
recovered for St 0. 

In some sense, the extra time derivative in the non- 
Fickian formulation plays a role similar to that of the dis- 
placement current in electrodynamics. In simulations of hy- 
dromagnetic fiows at low density, where the Alfven speed 
can be very large, the displacement current can be included 
with an artificially reduced value of the speed of light in 
order t o limit the Alfven speed to numerically acceptable 
values jMiller fc Stone 2600I) . 

A comment regarding centred finite difference schemes 
is here in order. In the steady state, the discretization of 
the cosmic ray diffusion model given by Eqs © and (|TTJ 
corresponds essentially to a conservative formulation of the 
diffusion term. (A conservative formulation involving a di- 
rect discretization of is not possible with a non-staggered 
mesh, because two first-order derivatives occur in two sep- 
arate equations.) As is well known, the discretization of 
the diffusion term on a centred non-staggered mesh means 
that structures at the mesh scale cannot be diffused (the 
discretization error for first derivatives becomes infinite). 
Therefore, we need to include weak Fickian diffusion in the 
cosmic ray energy equation. We refer to the corresponding 
(isotropic) diffusion coefficient as -/fpick, and it will be chosen 
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to be comparable to or less than the viscous and magnetic 
diffusivities. 

In the following we use the Pencil Code,^ a non- 
conservative, high-order, finite-difference code (sixth order 
in space and third order in time) for solving the compress- 
ible hydromagnetic equations. The non-Fickian diffusion for- 
mulation is invoked by using the cosmicrayf lux module. 
Whenever possible we display the results in non-dimensional 
form, normalizing in terms of physically relevant quantities. 
In all other cases we display the results in code units, which 
means that velocities are given in units of the sound speed 
Cs, length is given in units of k^^ (related to the scale of the 
box), density is given in units of the average density po, and 
magnetic field is given in units of ^fioPo Cs. The units of all 
other quantities can be worked out from these. For example, 
the unit of Qc is pocffci. For the interstellar medium with 
Po = 10"^" gcm"^, Cs = 10kms"\ and fei = 27r/100pc, the 
unit of the cosmic ray injection rate is 3 x lO"^'^ erg cm~^ s~^ , 
which is about 10% of the rate of energy injection by super- 
novae in the galactic disc (Mac Low & K lessen, ,200j). The 
unit for diffusivity is Csk^^ « 5 x 10^^ cm^ s^^. 

2.3 Cosmic ray diffusion near a magnetic X-point 

We test the field-aligned diffusion procedure by simulat- 
ing in two dimensions a magnetic field configuration similar 
to the X-point discussed in Sect. 12.21 In order to be able 
to impose normal-field boundary conditions, fi x B = 
at the domain boundaries, we modify the field to i? = 
(sin feix, — sin feij/, 0)"^, where fci is the smallest wavenum- 
ber in a periodic domain. So, for fci = 1 we consider the do- 
main — TT < {x, y) < n. The initial distribution of the cosmic 
ray energy density is ec — x, which has a constant gradient 
and therefore, with Fickian diffusion, Dc = V ■ (BB) would 
have a singularity initially. However, in the non-Fickian ap- 
proach Dc is not calculated as in Eq. JUJ, which resolves this 
problem. The evolution of ec for r = 0.1 is shown in Fig. |5| 
together with vectors showing the magnetic field. Note that 
the gradient of ec becomes small in the neighbourhood of the 
singularity of V ■ (BB) at the origin, so the otherwise singu- 
lar term that multiplies Vec has no effect on ec, as desired. 
In the case of the Fickian diffusion, the same final solution 
would have been obtained, but the initial reduction of the 
gradient in ec would have involved an infinitely large advec- 
tion speed Uc- In the non-Fickian approach, the maximum 

~ 1/2 

propagation speed is K^^ , thereby alleviating the numerical 
time step problem. 

Another example of field-aligned diffusion is shown in 
Fig. 13 where the magnetic field is given by S = Bq + 'V x A 
with Bo ~ O.lx and A = O.l£cos{kxx) cos{kyy) with 
kx — 4fci and ky — k\. Again, this magnetic field is held 
constant in time. The initial profile of ec oc exp(— r'^/2(T^), 
with = + (y-l-0.5)^, is a two-dimensional Gaussian of a 
half- width of cr = 0.07, positioned at (0, —0.5). We confirm 
that our implementation of cosmic ray diffusion allows us to 
model reliably rather complicated magnetic configurations. 
The lower panel of Fig. |5] confirms that, for large values of 
the Strouhal number, the wave nature of the telegraph equa- 
tion manifests itself and ec develops two waves propagating 

^ |http : //www. nordlta.dk/software/pencll- code I 
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Figure 2. Evolution of the cosmic ray energy density near a 
magnetic X-point: snapshots of ec (shown as contours and shades 
of grey/colour) for field-aligned diffusion along a fixed magnetic 
field B = (sinfcix, — sin fci j/, 0)"^ (shown as vectors) displayed for 
three times indicated at the top of each frame. 
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Figure 3. Magnetic field vectors togettier witli a grey/colour 

scale representation of ec in a kinematic calculation with 128^ 

1/2 

mesh points, for different values of t/ct, with K± = 0, K^^ = 

10^^, and ii'pick = 10"'^, at time t/r = 1 for two different values 
of T (=1 and 3, respectively). (Only part of the computational 
domain in the y direction is shown.) 



away from the initial maximum (similar to the dashed line 
in Fig. Pi. 



3 MACROSCOPIC EVOLUTION OF THE 
COSMIC RAY GAS 

3.1 Energy considerations 

In a closed domain, mass is conserved, i.e. (p) = po = 1, 
where angular brackets denote volume averaging. Then the 
hydromagnetic equations coupled with cosmic ray dynamics 
lead to the following set of equations for the cosmic ray 
energy Ec = (fie), the gas energy Eg = (eg), the kinetic 
energy i5k = {\pu^), and magnetic energy E^ ~ {B^) /2/io, 

^ = -Wc + {Qc), (15) 
dt 

dEn 



= Wc+Ws + Wn. + W{- (Qk>, 



dt 



(Qm). 



(16) 
(17) 
(18) 



Here, all the energies are referred to the unit volume. The 
terms Wc = (pcV ■u),Wg = (pgV ■ u), Wm = {u- {J x B)), 
and W{ = {u ■ f) result from work done against cosmic ray 
pressure, gas pressure, the Lorentz force, and the external 
forcing, respectively. Terms responsible for viscous and Joule 
heating and the cosmic ray energy source are simply given 
by the volume integrated terms in the original equations. 



Equations I15II - I18II imply that the total energy, iJtot = Ec + 
Eg + E]^ + Em, satisfies the simple conservation law 

d^tot 



dt 



(Qc) + W,. 



(19) 



Thus, the only sources of energy are the injection of cosmic 
rays and the external forcing of the turbulence. In the fol- 
lowing section we demonstrate how Ec can be enhanced by 
the conversion of kinetic energy. 



3.2 Compressional enhancement of cosmic ray 
energy 

We assume Qc = Wf = and that there is initially ki- 
netic energy that is later redistributed among gas and cosmic 
rays. We investigate, using a simple one-dimensional model 
{d/dy — d/dz — 0), how much energy can be converted 
into cosmic ray energy via the Wc term responsible for work 
done against cosmic ray pressure. As the initial condition, 
we use a sinusoidal perturbation of Ux and In p with unit 
amplitude and Ec = Eco = 1, £g = 1.8, and Ey, = 0.21. The 
evolution of velocity, cosmic ray and gas energies, as well as 
the entropy of the gas are shown in Fig. 0] Here the entropy 
s is defined as s = c„ ln(Cs /p^""*^), where = 7(7 — l)eg is 
the gas sound speed squared. It turns out that in this case 
about 78% of the kinetic energy is transformed into cosmic 
ray energy and only 22% into thermal energy. This result 
is however sensitive to the phase shift between density and 
velocity: if the density is initially uniform (keeping all other 
parameters unchanged) , the fractional energy going into cos- 
mic rays is only 23% and 77% go into thermal energy. 

These results demonstrate that, at least in principle, a 
sizable fraction of the kinetic energy can be converted into 
cosmic ray energy. Similar experiments have been made in 
earlier work with a similar model in the context of shock 
acceleration of cosmic rays ( see. e.g..lDrurY fc VollJll98ll: 
lJun et aPllQQ^) . In particular iKang fc Jone si Jl99Ci) showed 
that the efficiency of conversion varies strongly with 7c. 
However, the conversion of kinetic energy into cosmic ray 
energy requires a background of cosmic ray energy. De- 
creasing Ec from 1 to 0.1 lowers the fraction of compres- 
sionally produced cosmic ray energy density from 78% to 
21%. In contrast to dynamo theory where a weak seed mag- 
netic field is sufficient to produce equipartition magnetic 
fields (albeit only in three dimensions), there is no such 
mechanism for the cosmic ray energy. Thi s is related t o 
the anti-dynamo theorem for scalar fields (lKrauselll972l) . 
However, for three-dimensional compressible flows an ex- 
ponential dynamo-like amplification of a passive scalar is 
in principle possible if the passive scal ar is represented 
by inertial particles llElperin et alj[l996l) . Such a mecha- 
nism can work because inertial particles do not feel a pres- 
sure gradient. This c an lead to particle accumulation in 
temperature minima jElperin et alJll997^ and in vortices 
iBarge fc Somm erial |l995r ' Eod^onfc^^andenbur3 Il998l : 
Ijohansen et al.ll2004^ . However, in this paper cosmic ray 
particles are treated as non-inertial particles. 



3.3 Effect of cosmic ray pressure 

Cosmic rays can be confined at large scales by magnetic 
tension, where a strong magnetic field can more easily with- 
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Figure 4. Velocity, cosmic ray and gas energy densities, and 
entropy in an experiment with a non-linear sound wave that piles 
up to a shock (7c = 5/3) . Note the significant conversion of kinetic 
energy into cosmic ray energy. The conversion into gas energy 
is comparatively small even though there is noticeable entropy 
enhancement due to the shock. Curves obtained for different times 
are shown with different lino types as labelled in the first panel. 
Time is given in units of fc;["''"(i?co/po)~"^''^- 

stand deformation driven by cosmic ray pressure gradi- 
ents. This could provide a natural mechanism for producing 
equipartition between cosmic rays and the magnetic field. 
This feature can be simulated in two dimensions in a dou- 
bly periodic domain — tt < (x,y) < n, with ki — 1. The 
results are illustrated in Fig. |^ where we have a magnetic 
tube in 1 < 3/ < 2 with its axis along the x direction. We 
have implemented two local cosmic ray sources with energy 
injection profiles 
2 

Oc = QcoX]exp{-5 [x^ + {y-y^f]/R^}, (20) 
1=1 

i.e., both located on the y axis, centred at j/i =0 and 
J/2 = 7r/2; the initial half-widths for both sources is R — 0.13, 
so that one source is within the magnetic tube and the other, 
outside it. In this experiment, cosmic ray diffusion is negli- 
gible (A'li — K± — and ifpick ~ 0.01) as we intend to ex- 
plore the effects of cosmic ray pressure alone. As expected, 
expansion proceeds nearly isotropically outside the magnetic 
structure, but the cosmic ray energy density is channelled 
preferentially along field lines inside the tube. At the end 
of the run, the aspect ratio of the cosmic ray distribution 
is about two to one inside the tube. For values of Qc sig- 
nificantly larger than about 10, the gas density decreases 



t Cs ki = 0.5 




0.0 1.0 2.0 3.0 4.0 

« Cs fci = 2.0 




t Cs ki = 10.0 




Figure 5. Cosmic ray energy density at times indicated at the 
top of each panel. Cosmic rays expand from two sources (with 
injection rate Qc = 10 for each), one inside a magnetic flux tube 
and the other one outside. Magnetic lines are shown with white 
solid curves whose density is proportional to the field strength. 

strongly so as to maintain pressure equilibrium and oppose 
expansion driven by cosmic rays. 

This confirms that cosmic ray dynamics can be strongly 
affected by the approximate pressure balance in the ISM. 

3.4 Cosmic rays in a partially ordered magnetic 
field 

In this section we briefiy explore the effects of a random mag- 
netic field on the evolution of the cosmic ray gas. A random 
component of the interstellar magnetic field can facilitate 
the isotropic spreading of cosmic rays across the large-scale, 
preferentially horizontal magnetic field in the Galactic disc. 
In addition, a turbulent magnetic field can enhance cosmic 
ray diffusion by destroying the com pound diffusion effect 
l|Ptuskinll979HK6ta fc Jokioi i 2000, and references therein) 
due to the exponential local divergence of magnetic lines. 

To allow for cosmic ray losses through the x bound- 
aries, we relax the assumption of periodicity in that di- 
rection. At a; = ±7r, we assume Cc = 0, together with 
dp/dx — de^/dx — 0. This implies that cosmic rays may 
be lost from the domain but gas may not. In the y direction 
we again use periodic boundary conditions. 

We consider a two-dimensional system with a regular 
magnetic field Bg directed along the j/-axis and confined to 
a flux tube as shown in Fig.|Hl where the field strength has a 
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Figure 6. Cosmic ray energy density (colour/grey scale coded, 
witli darker/blue sliades corresponding to smaller values) together 
with magnetic field lines (soUd) in a two-dimensional simulation 
with a fixed magnetic flux tube centred around x = 1.5 and a 
random magnetic field superimposed on it. Here, K^^ = 0.1, = 
0, and T = 3. 




Figure 7. Cosmic ray energy density from the model of Fig. ISl 
averaged in the y direction for times 125 X 2" with n = 0, 4. 
The magnetic tube is located at x = 1.5 leading to an asymmetric 
distribution of cosmic ray energy density. 



profile Bo oc sech'^[(a;— 1.5)/0.5]. An isotropic random mag- 
netic field SB is superimposed on Bo, with SB'^/Bl = 1 at 
a; = 1.5 where Bq is maximum; the magnetic field does not 
evolve. The random magnetic field is implemented in terms 
of a magnetic vector potential given as white noise with 
Gaussian probability density which, because of two dimen- 
sions, implies a power spectrum for the magnetic energy. 
We also assume zero velocity for all times, so we just advance 
Eqs and (IIH in time, using Eq. JSJ. Cosmic rays are in- 
jected at a constant rate across the domain, Qc = const. 

In Fig. il we show the result of such a calculation with 
K± = ; the distribution of cosmic rays in x is asymmetric 
reflecting the asymmetry in the relative amount of disorder 
of the magnetic field, SB^/B^. This asymmetry can be seen 
more clearly in Fig. Q which shows the evolution of cosmic 
ray energy density averaged in the j/-direction. (Note how- 




Figure 8. The profile of K^f^'' (solid) obtained from Eq. 1221 
using Ec corresponding to the upper curve of Fig. |7] and x = 
3B^/B^ (dashed), where By has both large-scale and random 
parts, whereas B^ is a purely random magnetic field. Here, 



■ 0.023 is the maximum value of K 



(eff) 



ever that the steady state is only attained after very long 
times. Here, t = 2000 corresponds to trK^^kf — 600.) 

The effective perpendicular diffusivity due to the ran- 
domness of the magnetic field, K'f^\x), can be obtained 
from the steady-state equation 



(21) 



which can be integrated to obtain 



(22) 



where xq is the position where dec/dx — 0. The resulting 



(off) 



profile of K 



shown in Fig. |H| along with 



(23) 



confirms that the effective perpendicular diffusion is con- 
trolled b y the degree of randomness of the magnetic field 
(see, e.g., IChuvikin fc Ptuskinliriol . 



4 COSMIC RAYS IN A MAGNETIC FIELD 
PRODUCED BY DYNAMO ACTION 

Three-dimensional turbulence is capable of dynamo action 
for sufficiently large magnetic Reynolds numbers, and 
the dynamo-generated magnetic fiel d organizes itse l f into 
random flux tubes or sh eet s fe.g. |ZeMo^c^^^lJ ^90|: 
i Braridenburg et all 119951 : iBrandenburg fc SubramaniarJ 
I2OO5I and references therein). Most studies of cosmic ray 
dynamics neglect the speciflc features of the magnetic 
flelds produced by turbulent dynamos. We provide here a 
preliminary discussion of cosmic ray evolution in a magnetic 
held generated by a turbulent flow of electrically conducting 
fluid. The magnetic held structure of these simulations is 
realistic enough to include important physical effects such 
as the enhancement of cosmic ray diffusion by turbulent 
flelds, as mentioned in Sect. 13.41 

Magnetic field produced by the dynamo action is rather 
different from that prescribed as, say, a random vector field 
with given spectrum and Gaussian statistical properties of 
the components. In contrast to such ad hoc models, dynamo 
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magnetic fields can be strongly intermittent (i.e., dominated 
by intense magne tic filaments, ribbons and sheets) and vary- 
ing in time (see iBrandenburg fc SubramanianI l2005l and 
references therein); both features can affect the propaga- 
tion of charged particles. Moreover, since both gas fiow and 
magnetic field are random (in space and time), any relation 
between cosmic ray energy density and other parameters of 
the medium (e.g., magnetic energy density or gas density) 
can only be statistical. Therefore, we expect that the en- 
ergy density of cosmic rays can locally (and at any given 
moment) significantly exceed, say, the magnetic energy den- 
sity. However, one would expect that some form of equiparti- 
tion between energy densities (or forces due to) cosmic rays 
and magnetic fields can be maintained on average. We note, 
however, that simul ations have not fully co nfirmed these ex- 
pectations; see also lPadoan fc Scald ll2005l) . 

Our model is realistic with respect to modelling fully 
nonlinear dynamo action as we simulate consistently both 
a randomly forced fiow and the magnetic field produced by 
it, by solving both the Navier-Stokes and induction equa- 
tions (with the Lorentz force included in the former, and the 
velocity field obtained from the Navier-Stokes equation in 
the latter). The turbulent motions in our model axe driven 
by a random force explicitly included in the Navier-Stokes 
equation. In reality, interstellar turbulence is driven by su- 
pernova explosions that produce strongly compressible flows 
with very large local Mach numbers l ocally (some aspects 
of th e relevant models are reviewed bv lMac Low fc KlessenI 
|2004) . However, we deliberately restrain ourselves from a de- 
tailed discussion of such more realistic models here (which 
would also include stratification, disc-halo connections, ve- 
locity shear, etc.), but instead explore just the effects of 
magnetic intermittency and variability. We believe that our 
simulations capture at least some of the most important 
effects of interstellar dynamo action on the cosmic ray prop- 
agation (within the limits of our model for the cosmic rays). 

The turbulence in our simulations is driven helically by 
a forcing function / in t he Navier-Stokes equ ation, as was 
done in the simulations of lBrandenburd (|2Q^), for example. 
At a; = ±7r, we use stres s-free normal field boundary c ondi- 
tions (as was also done in lBrandenbxirg fc Doblerl2001^ . and 
assume Cc = on the boundaries as in Sect. 13.41 In the other 
directions we take periodic boundary conditions. Our anal- 
ysis of the results presented below only uses positions that 
are some distance away from the domain boundaries {L^/S 
on both boundaries) to reduce their infiuence. (Including 
boundary points merely tends to decrease the magnitude of 
the correlation coefficients between the various energy den- 
sities, but it does not change the results qualitatively.) The 
forcing function is given in Appendix |B| and its (dimension- 
less) amplitude for the simulation shown here is chosen to 
be fo — 2, which produces an rms Mach number of about 
1.2. 

The forcing wavenumber is chosen to be fcf — 1.5 fci. 
This value is close to the wavenumber corresponding to the 
box size, ki = 271/ Lx, so we do not expect to have clearly dis- 
tinct large-scale and small-scale magnetic fields. Generally, 
the flow helicity allows us to obtain dynamo action at rela- 
tively small values of the magnetic Reynolds number deflned 
as -Rm = Urins/{rjk{). However, because of the non-periodic 
boundaries in the x direction, and also because of the weak 
scale separation (fcf/fci is not very large), the critical value 




Figure 9. Time series of magnetic (i?m), kinetic (i?k) '^'nd cosmic 
ray (-Ec) energies in a dynamo simulation. Here, time is given 
in turnover times (urms^f)"^, and ScO = L^Qc/K^^ is used to 
normalize energies per unit volume. The thermal energy of the 
gas is constant with Eg/ccO ~ 0.7. 



of -Rm with respect to the onset of dynamo action is still 
around i?m,cr = 30, which is similar to what would be ex- 
pected for a non-helical random flow in a periodic domain. 
The simulation presented here has Rm ~ 150. The kinematic 
growth rate of the rms magnetic fleld is about 0.06 Urmsfcf . 
In Fig. 1^ we show the evolution of the magnetic energy to- 
gether with kinetic and cosmic ray energies. We see that 
the magnetic fleld grows exponentially for t ^ 150/(urmsfcf) 
and then saturates — in agreement with earlier simulations 
quoted above. We note that the energy density of cosmic 
rays is much larger than magnetic energy density at these 
early times; nevertheless, the cosmic ray energy increases 
rather slowly after t <: 50/(urmsfcf). The steady-state energy 
density of cosmic rays is controlled by their injection rate Qc 
and their diffusivity: solutions of Eq. 12H are proportional 
to Qc/K'f'^^ ■ However, the effective diffusivity of cosmic rays 
is controlled by the degree of tangling of the magnetic fleld 
rather than by the fleld strength itself; see, e.g., Eq. (1231 . It 
is not surprising, then, that even a weak magnetic fleld can 
conflne cosmic rays at early times in this model. The lin- 
ear dependence of the steady-state energy density of cosmic 
rays on their injection rate is a direct consequence of the (al- 
most) linear nature of the cosmic ray dynamics as described 
by Eq. (|5J ; the only nonlinearity here is that the cosmic ray 
energy density affects the flow through the pressure term, 
and then the velocity field enters the induction equation and 
the advection term for the cosmic rays. However, this non- 
linearity is not very strong, and our simulations confirm a 
linear dependence of ec on Qc within a broad range of the 
latter (at least two orders of magnitude) . The magnetic field 
part Bo is understood, in the present context, as an average 
over a scale smaller than the domain size but larger than, 
say, the gyro radius of cosmic ray particles. 

For the simulation shown here we have chosen Qc = 
0.01, which yields a steady-state cosmic ray energy of -Ec ~ 1 
in units of L^Qc/Kp The other parameters of the simula- 
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Figure 10. Cosmic ray energy density (colour/grcy-scale coded, 
with lighter shades/redder colour corresponding to larger values) 
and magnetic field vectors in a slice taken from a dynamo simula- 
tion. The magnetic field vectors are more numerous where mag- 
netic field strength is larger. 



tion presented here are K± = 0, K^^ = 0.3, Kpick ~ 2x 10~^, 
T = 0.3, = 5 X 10"'^, u = 0.5. Furthermore, because the 
Mach number is slightly larger than unity, an additional bulk 
viscosity proportional to the negative velocity divergence has 
been i ncluded. This is usua lly referred to as a shock viscos- 
ity; see lHaueen et all (|200^ for details and the definition of 
a non-dimensional parameter Cshock which is here chosen to 
be 10. The value of A'n is chosen to be close to the maxi- 
mum Alfven speed squared. The magnetic field produced by 
the dynamo has pronounced magnetic filaments whose half- 
width (radius) is about £ — 0.2, which is consistent with the 
estimate £ ~ Tvk^^ Rm[c^ suggested bv nSubramanian. 
For T = 0.3 and £ = 0.2, we have St « 1 from Eq. JTSJ. The 
steady-state mean kinetic energy density depends directly 
on the intensity of the forcing. On the other hand, the ra- 
tio of magnetic to kinetic energy densities is controlled by 
the nature of the dynamo action. The above parameter val- 
ues have been chosen as to ensure that the energy densities 
of magnetic field and cosmic rays are of the same order of 
magnitude in the statistically steady state. 

In Fig. 1101 we show a typical cross-section of the cos- 
mic ray energy density and magnetic field vectors from 
the three-dimensional dynamo simulation of Fig. |5| at t = 
250/('Urmsfcf). The cosmic ray energy density declines to- 
ward the boundaries at a; = ±n, where the boundary con- 
dition Cc = is imposed, and shows some moderate vari- 
ation inside the domain. There is no pronounced correla- 
tion with magnetic field strength even though imprints of 
the field-aligned diffusion can clearly be seen, e.g., between 
{x,y)ki = (-1,-1) and (0,0). We show in Fig. [TTla two- 
dimensional joint probability density function of log and 



o 





-2 

-4 

-6 

-8 
-10 




5f i» 



0.2 0.4 0.6 0.8 1.0 1.2 1.4 

Sc / EcO 



0.0 0.0003 0.009 0.06 0.3 0.8 
p(ec, logB^) 

Figure 11. Two-dimensional histogram (or joint probability 
density) of magnetic pressure and cosmic ray energy density. Here, 
Bco = l^xQc/^W is used to normalize Sc- The two-dimensional 
probability density is calculated using only points at a distance 
greater than L^/S from the boundaries in an attempt to avoid the 
regions where the distribution of ec is affected by the boundary 
conditions. 
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Figure 12. As in Fig. 1111 but for gas density and cosmic ray 
energy density, showing a modest correlation between the two. 
The correlation coefficient is r = 0.54, and the straight line is a 
best-fitting line. 



Be (normalized to unit integral as usual), which demon- 
strates the lack of any noticeable correlation between these 
variables. The finite lifetime of magnetic structures pro- 
duced by the dynamo must be one of the reasons of the 
lack of correlation between the two variables. There is some 
correlation between gas density and cosmic ray energy den- 
sity, as shown in Fig. 1121 but the cross-correlation coefficient 
is only 0.54, with the best-fitting dependence Cc/eco — p/po- 
If the injection rate of cosmic rays is reduced by a factor 
of ten to Qc = 10"'^, the resulting steady-state mean value 
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of the cosmic ray energy density is found to be reduced by 
about the same factor. The relation between cosmic ray en- 
ergy density and gas density still appears to be nearly linear, 
but the cross-correlation coefficient is now larger, varying 
with time in the range 0.7-1. 

We have also explored a run with a Mach number of 
0.2 (achieved by using a weaker driving force; /o = 0.05) 
but with the same injection rate of cosmic rays as above, 
Qc = 0.01. The resulting steady-state mean energy density 
of cosmic rays exceeds those of magnetic field and turbu- 
lence, Sc/eco 1, -Ek/eco ==; 0.02 and Sm/eco ==; 0.01. This 
produces significant anticorrelation between cosmic ray en- 
ergy density and gas density (cross-correlation coefficient of 
—0.94), with a linear dependence between Cc and p. 

The latter anticorrelation may be attributed to the av- 
erage pressure equilibrium in the domain, while a positive 
correlation in the supersonic flow may arise as both gas and 
cosmic rays are compressed by the gas flow. We have con- 
flrmed that no positive correlation between cosmic rays and 
gas density occurs if the cosmic ray advection is neglected. 

The model illustrated in Figs |U] and IIUI is close to en- 
ergy equipartition between cosmic rays, magnetic fleld and 
turbulence. We note however that the Lorentz force and 
the cosmic ray pressure gradient have very different magni- 
tudes because the field-aligned cosmic ray diffusivity is much 
larger than the magnetic diffusivity. As a result, the cosmic 
rays are distributed more uniformly than the magnetic fleld 
and the gas density and so the cosmic ray pressure gradi- 
ent is comparatively small. For the values of the diffusivities 
given above, the ratio the rms cosmic ray pressure gradi- 
ent, Fc, and the rms Lorentz force Fm, is typically about 
0.1 of the ratio of the corresponding mean energy densities; 
this also applies if the Lorentz force is replaced by the gra- 
dient of kinetic energy density. The typical length scale of 
the magnetic field is about I ~ loR^^l^ and F^ ~ em/^, 
with Zo — 100 pc the turbulent scale and -Rm.cr ~ 30 the 
critical magnetic Reynolds number for the onset of dynamo 
action (see above). The length scale of the cosmic ray dis- 
tribution can be estimated as the diffusion scale over the 
confinement time Tc ~ 10^ yr, Ic — {K^TcY^'^ . Then, for 
7^11 = lO^^cm^s-i, 

Fm Ic Em 30 -Bm 

This conclusion appears to be model-independent and sug- 
gests that energy equipartition between cosmic rays and 
other constituents of the interstellar medium does not nec- 
essarily imply that cosmic rays play an important role in the 
dynamical balance. 



5 DISCUSSION AND CONCLUSIONS 

We have presented a preliminary analysis of cosmic ray prop- 
agation in a magnetic field produced by dynamo action of 
a turbulent flow. The conflnement of cosmic rays resulting 
from their scattering by magnetohydrodynamic waves can 
be modelled with an equation similar to Eq. 0, where the 
advection veloci ty is a linear c ombination of gas velocity and 
Alfven velocity JSkillindll975l) . Our results are based on ad- 
vection with the local gas velocity, ^adoan & Scalo (200^ 
considered local variations in cosmic ray density in the case 



where the advection velocity is given by the Alfven veloc- 

1/2 

ity. They predict that ec oc , with m the ion density. 
This scaling is expected if the diffusive streaming velocity, 
— 7^11 Vec, and the effects of cosmic ray pressure are negligi- 
ble. Our model can be adapted to test and generalize these 
results; the anticorrelation between Cc and gas density in one 
of our models (with low Mach number) seems to be a direct 
consequence of pressure balance, while a positive correla- 
tion (obtained at larger Mach number) may reflect the fact 
that both cosmic rays and thermal gas experience similar 
compression by the gas fiow. We have shown that our model 
captures naturally the dependence of the effective diffusivity 
of cosmic rays on the ratio of random to ordered magnetic 
fleld, 5B^/Bi. 

The diffusivity of cosmic rays along the magnetic fleld is 
rather large; the corresponding Strouhal number, defined in 
Eq. 1131 may significantly exceed unity, as shown in Eq. 11411 . 
For comparison, a similar estimate yields St ~ 1 for the tur- 
bulent kinetic and magnetic diffusivities in the ISM. This 
motivates our suggestion that the standard Fickian diffu- 
sion model, which leads to the classical diffusion equation, 
may be a poor approximation for cosmic rays, and a more 
accurate description leading to some form of the telegraph 
equation might be more appropriate. Formally, the differ- 
ence between the two approximations consists of retaining, 
in the latter approximation, higher-order terms in the corre- 
lation time of the random process underlying diffusion. We 
have introduced this effect to alleviate numerical problems, 
but it can be a real physical effect which deserves further 
careful study. 

In summary, we have found that the cosmic ray dis- 
tribution can be more uniform than the distributions of 
magnetic field and gas density. Consequently, we may ar- 
gue that energy equipartition between cosmic rays and other 
constituents of the interstellar medium does not necessarily 
imply that cosmic rays play a significant role in the dynam- 
ical balance. 



APPENDIX A: BOUNDEDNESS OF COSMIC 
RAY ENERGY DENSITY 

In this section we show that, in a closed or periodic domain, 
max(ec) can only decrease as a result of (tensorial) diffusion. 
This is useful for showing that the diverging behaviour of Uc 
does not produce a singularity in Cc', cf. Sect. 12.21 In order 
to avoid interference from other effects, we assume that the 
evolution of ec is only governed by diffusion, i.e. 

^=V,(7^,,V,ec). (Al) 

Note also that max(ec) — (e")^'^" for n oo. Here, angular 
brackets denote volume averages. Thus, using integration by 
parts, we have 

= ^n{n - 1) (er'i^»,(V.ec)(V,ec)) 

< (for any value of n > 1). (A2) 

The last inequality assumes that the diffusion tensor is pos- 
itive definite, which is true in our case, because 

K,j{\7,e,){\7Je,) = {K^^-K^){B ■Ve.f+K^iVe.f (A3) 
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is positive. Therefore, max(ec) must decrease with time. 



APPENDIX B: THE FORCING FUNCTION 

For completeness we specify here the forcing function used 
in the present paper^. It is defined as 



fix, t) = Re{Nfk(t) exp[ifc(f) • x + icj){t)]}, 



(Bl) 



where x is the position vector. The wavevector k{t) and the 
random phase — tt < 4'{t) ^ tt change at every time step, 
so f{x,t) is J-correlated in time. For the time-integrated 
forcing function to be independent of the length of the 
time step 5t, the normalization factor N has to be propor- 
tional to St~^^^ . On dimensional grounds it is chosen to be 
N = fopoCs{\k\cs/Sty^\ where /o is a dimensionless forc- 
ing amplitude. At each time step we select randomly one of 
many possible wavevectors in a certain range around a given 
forcing wavenumber. The average wavenumber is referred to 
as k{. Two different wavenumber intervals are considered: 
1-2 for k{ = 1.5 and 4.5-5.5 for k{ = 5. We force the system 
with transverse helical waves, 



/fe = R ■ j'""*^"'-* 



.,1 ri ^^j ^^^ijkkk 

With Rij = 



where a = 1 for positive helicity of the forcing function 



(B2) 



(B3) 



is a non-helical forcing function, and e is an arbitrary unit 
vector not aligned with fe; note that |/fe| =1. 
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